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Abstract 



We study the motion of light in the gravitational field of two Schwarzschild 
black holes, making the approximation that they are far apart, so that the 
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i^ ' motion of light rays in the neighborhood of one black hole can be consid- 



^^ ' ered to be the result of the action of each black hole separately. Using this 

^ ■ approximation, the dynamics is reduced to a 2-dimensional map, which we 

^ ' study both numerically and analytically. The map is found to be chaotic, 

with a fractal basin boundary separating the possible outcomes of the orbits 
(escape or falling into one of the black holes). In the limit of large separation 
distances, the basin boundary becomes a self-similar Cantor set, and we find 
that the box-counting dimension decays slowly with the separation distance, 
following a logarithmic decay law. 
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I. INTRODUCTION 

In this article, we study the motion of hght (null geodesies) in the gravitational field of two 
non-rotating Schwarzschild black holes. In general relativity, solutions of the field equations 
describing more than one purely gravitational sources are necessarily non- stationary, because 
gravity is always attractive (we are not considering exotic matter); there is no possibility of 
arbitrarily "pinning" sources as is done in Newtonian gravitation, because of the automatic 
self-consistency of the nonlinear Einstein's equations. If we demand that the two black 
holes be fixed in space, then the solution includes a conical singularity (a "strut") lying 
on the axis on which the two masses are located fl]]. This singularity appears as a natural 
consequence of the field equations, and it is necessary to keep the two masses from falling 
towards each other. However, this singularity would have to be made of very exotic matter, 
and this solution does not describe any realistic system in astrophysics. The real solution 
to the relativistic two-body problem can have no conical singularities, and it is necessarily 
non-stationary. The two black holes will spiral around each other, emitting gravitational 
waves, which makes this problem even more difficult. There is no exact solution for the 
relativistic two-body problem, and even a numerical solution has eluded the most powerful 
computers. 

In order to cope with this problem, Contopoulos and others 0,^] have used the 
Majumdar-Papetrou solution |l[| to study the dynamics of test particles in a space-time 
with two black holes. The Majumdar-Papetrou metric used by Contopoulos describes two 
non- rotating black holes with extreme electric charge {Q = M in relativistic units), whose 
gravitational pull is exactly matched by their electrostatic repulsion, thereby allowing a 
static mass configuration. They have found that in this metric the motion of both light and 
massive particles is chaotic, with a fractal invariant set and a fractal basin boundary. How- 
ever, it is very unlikely that the Majumdar-Papetrou metric describes realistic astronomical 
objects, since there is no known realistic astrophysical process by which a black hole with 
extreme charge could be formed. Even though the two black holes with extreme charge have 



proven useful as toy models, it is important to address the more realistic problem of two 
uncharged black holes. This is what we do in this article, for the motion of light and other 
massless particles. 

In order to overcome the fact that there is no static solution for the two-black-hole system, 
we consider the case when the two black holes are far apart, with a distance much larger than 
their Schwarzschild radii. In this case, nonlinear effects in the field equations are expected 
to be small, and we can approximate the motion of test particles in the neighborhood of 
one of the masses as being the result of the field of that mass alone, and disregard the 
effect of the other black hole as being negligible. Using this approximation, the motion of 
test particles in the two-black-hole system is treated as a combination of motions caused 
by isolated Schwarzschild black holes. Since the equations of motion for the Schwarzschild 
geometry can be analytically integrated, our dynamical system is reduced to a map, which is 
much easier to study than a system of ordinary differential equations. This scattering map 
is built in Section 2, for the simple case of two black holes with equal masses. In Section 3, 
we show that this map has a fractal basin boundary separating the possible outcomes of a 
light ray in the two-black-hole field, namely, falling into either of the black holes or escaping 
towards the asymptotically plane infinity. The fractal (box-counting) dimension of this basin 
boundary is numerically calculated, and the sensitivity to initial conditions implied by the 
fractal nature of the boundary is thereby quantified. In Section 4 we use explicitly the 
condition of large separation between the black holes. In this limit, the basin boundary 
becomes a self-similar Cantor set, which allows us to obtain some analytical results. One 
of our main results is that the fractal dimension decays very slowly (logarithmically) with 
distance. The slow decay of the fractal dimension makes it more likely that the fractal nature 
of the basin boundary has some importance for astrophysics. In section 5, we consider the 
case of two black holes with unequal masses, in the limit of a large separation; we find that 
the logarithmic decay law of the fractal dimension for large distances is also valid in this 
case. In section 6, we summarize our results and draw some conclusions. 



II. SCATTERING MAP FOR TWO BLACK HOLES WITH EQUAL MASSES 

We begin by reviewing some basic results concerning the motion of test particles in the 
field of an isolated Schwarzschild black hole [|,Q . We consider specifically the case of null 
geodesies, which concerns us most, but many features of the dynamics also apply to massive 
particles. 

The Schwarzschild metric is written in spherical coordinates as: 

ds^ = f 1 jdf ^-r^dQ^, 

r 

with dfi^ = dO^ + sin^ 6'(i0^ being the element of unit area, and t is the time measured from 
a distant observer. M is the black hole's mass in geometrized units. We are interested only 
in the region of space-time outside the event horizon, r > 2M. Due to the conservation 
of angular momentum, test particles move on a plane, which can be chosen as ^ = 7r/2. 
The plane whereon the motion occurs is then described by the coordinates r and 0. The 
geodesic equations which describe trajectories of test particles on this plane can be analyti- 
cally integrated by means of elliptical functions 0] . Here we are interested in the scattering 
of null geodesies by the black hole. A light ray coming from infinity towards the black hole 
is characterized by the impact parameter h defined by the ratio h = L/E, where the angular 
momentum L and the energy E are constants of motion given by: 

E = 
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A is the geodesic's affine parameter. For null geodesies only the ratio of L and E is of 
importance to the dynamics. In the asymptotically plane region r — > oo, h corresponds to 
the usual impact parameter of classical scattering problems. 

If the impact parameter is below the critical value h^ = 3^/3M, the trajectory of the light 
ray spirals down the event horizon and plunges into the black hole. If 6 > be, the trajectory 



circles the black hole and escapes again towards infinity, being defiected by an angle A. The 
lowest value P of the radial coordinate r along the trajectory (the "perihelium" ) is given 
by: 

■^ (1) 



P-2M 
Following Chandrasekhar [Q, we define the quantities Q, k and x by: 

Q2 = (P-2M)(P + 6M); (2) 



2 Q-P + QM 



2Q 



(3) 



2/ /^N Q-P + 2M 



The scattering angle A is then given by 

A = 7r-2/(6), (5) 

and the function f{b) is: 

fib) = 2^[Kik)-F{x/2,k)]. (6) 

Here F is the Jacobian elliptic integral and K is the complete elliptic integral. In Fig. 
1 we show a plot of A (6). As b approaches the critical value b^ from above, A goes to 
infinity; trajectories with b sufficiently near b^ can circle the black hole an arbitrary number 
of times before escaping, and for b = be, the light ray makes an infinite number of rotations, 
and never escapes. This is a consequence of the existence of an unstable periodic orbit at 
r = 3M, which appears as a maximum in the effective potential. The orbits with b = be 
spiral towards the r = 3M orbit, and in the language of dynamical systems they make up 
the stable manifold associated with this periodic orbit. 

The fact that A assumes values above vr for a non-zero range of b implies the existence 
of a rainbow singularity in the scattering cross section; this is to be contrasted with the 



Newtonian Rutherford scattering, which shows no such singularities. In fact, A assumes 
arbitrarily large values, and the differential cross section at any given angle 9 is made up 
by an infinite number of contributions arising from trajectories with A = 6, A = 6 -\- 27r, in 
general A = ^ + 2mT, corresponding to trajectories that circle the black hole n times before 
being scattered towards 6. However, large values of n correspond to very low ranges of b: the 
set of trajectories that scatters hj 6 + 2mT has a measure that decreases very rapidly with 
n. Chandrasekhar [^ shows that the impact parameter bn corresponding to a scattering by 
6 + 2mT for large values of n is given approximately by: 

&„ = 6e + 3.48Me-(''+2n-)_ (7) 

This expression shows that the measure of the set of trajectories scattered by 6 + 2mT decays 
exponentially with n, and the contribution of orbits with large n to the cross section is small. 
In fact, we shall see later that in many cases it is a good approximation to consider only 
orbits with n = 0. 

After reviewing some properties of an isolated black hole, we now consider the case of 
two black holes with equal mass M (we consider the case of different masses in Section 5). 
As we mentioned in the introduction, there is no exact solution of Einstein's field equations 
that describes this system. Because of this, we assume that the two black holes are separated 
by a distance D much larger than their Schwarzschild radius 2M; in this limit the nonlinear 
interaction between the two gravitational fields can be ignored. In a real system, the two 
black holes will be rotating around their center of mass; however, their rotation speed is 
much smaller than the velocity of light. We can thus consider the two black holes to be 
fixed in space, without incurring in too much error. Notice that this approximation might 
not be valid for massive test particles. 

We are interested in the orbits that never escape to infinity nor fall into one of the event 
horizons; these orbits make up the basin boundary of the system, which will be discussed 
later in more detail. For the orbits not to escape, they need to have impact parameters such 
that they are scattered by at least vr by one of the black holes. In the case of an isolated 



black hole, this corresponds to an impact parameter lower than b = besc ~ 5.35696M, which 
is less that three times the Schwarzschild radius. Since in our approximation D ^ 2M, for 
the purpose of finding the basin boundary we can consider that the light rays are scattered 
by each black hole separately, the other black hole being too far away to make a significant 
difference in the scattering. After suffering a scattering by one of the black holes, the ray may 
reach the other black hole, depending on its emerging trajectory after the first scattering. It 
is then scattered again, and may return to the first black hole, and so on. Since D 3> 2M , 
we consider the scattering process of each black hole separately and use formulas (|^) and 
(P) to determine the deflection angle due to each black hole as a function of the incident 
impact parameter. 

By making the approximations mentioned above, we reduce the motion of light in the 
two-black-hole space-time to a 2-d map, as has been done in to study general features of 
chaotic scattering. To do this, we make the further assumption that the light rays have zero 
angular momentum in the direction of the axial symmetry axis, on which lie the two black 
holes; the orbits are then confined to a plane containing the two black holes. Due to the 
axial symmetry of the system, the motions on all such planes are similar. Now suppose we 
have a light ray escaping from one of the black holes with impact parameter 6„ and with an 
escaping angle 0„ with respect to the symmetry axis, as shown schematically in Figure 2. 
Since the black holes are considered to be very far apart, the impact parameter 6„+i of the 
light ray with respect to the other black hole is the segment / shown in Figure 2 (one black 
hole can be considered to be "at infinity" as regards the other). We use the convention that 
positive values of b means that the ray is directed to the right side of the black hole, and rays 
with negative b are directed to the left. From elementary geometry, we have I = bn + D sin 0. 
The deflection angle is given by A(6„+i). The map is then written as: 

bn+l = bn + DsilKpn, (8) 

4>n+i = TT + 0„ - A(6„+i) mod 27r. (9) 



The angles 0„ are measured counterclockwisely with respect to each black hole; the first 
term in Eq. (^ comes from the change in the angle's orientation necessary to take account 
of that. 

Consider the initial conditions 60 = ^esc and (/)o = 0. Since A{besc) = tt, we see from the 
above Equations that these values of b and are a fixed point of the map. It corresponds 
to the periodic orbit depicted in Figure 3a, which revolves around the black holes, making 
a U-turn at each black hole and then heading towards the other. Another periodic orbit is 
shown in Fig. 3b. This orbit is such that bn+i = —bn and (pn+i = —<Pn- Inserting these 
conditions in Eqs. (H) and @, we find 26o = — -Dsin^o and A(6o) = vr + 2\(J)q\ (remember 
that the angles are defined modulus 27r), with 69 > 0. 0o is given by the solution of the 
equation 

D 



A(^-sin|0o|) =7r + 2| 



These are the simplest periodic orbits, but there are many others. 

We observe that Equations (H) and @) are valid only as long as b remains within the 
range be < b < besc- If b falls out of this interval, the ray either escapes or falls into one of 
the black holes, and the iteration must be stopped. 

III. ANALYSIS OF THE SCATTERING MAP 

We now proceed to study in detail the map defined by Equations (|) and @. We begin 
by a direct numerical investigation of these equations. 

In order to iterate Equations (^ and (|^) for given initial values 0o and 60, we first have 
to be able to calculate the deflection angle A for a given impact parameter b. To do this, 
we must begin by finding the "perihelium distance" P corresponding to b; this is done by 
solving the third-order equation (|1]) for P. We use the well-known Newton-Raphson method, 
which guarantees a very fast convergence [^. We then use Eqs. (^, (D and (^) to calculate 
Q, k and X; and we finally substitute these quantities in Eqs. (^ and (^ to obtain A(6). 
The elliptical functions F and K are computed by numerical routines found in [0. 

8 



Depending on its initial conditions, a light ray may either fall into one black hole, fall 
into the other black hole, or escape towards infinity. The set of initial conditions which leads 
to each of these outcomes is called the basin of that outcome. In our numerical iteration of 
the map (§,^, we are interested in obtaining a basin portrait of the system. To do this, we 
have to choose a set of initial conditions and iterate them to find out to which basin they 
belong. Our choice is the one-dimensional set with 0o = and an interval of b. As we have 
seen in the previous Section, if \b\ < be, the light ray always falls into the event horizon, and 
if \b\ > besc, it always escapes. We thus choose the interval to he be < b < besc- We divide 
this segment into 5000 points, and iterate the map (H,P) for each of these initial conditions, 
recording the final outcome for each point: if at any point in the iteration |6„| < be, this 
means that the light ray falls into one of the black holes, and if |&„| > bese, it escapes to 
infinity. We define the discrete-valued function g{b) to be 1 if the orbit with initial conditions 
(po = 0, bo = b falls into one of the black holes, -1 if it falls into the other, and if it escapes 
to the asymptotically plane region. g{b) gives a picture of the intersections of the three 
basins with the segment 0o = 0, &c < ^ < bese- 

The result is shown in Fig. 4a, for D = 15M. We see that there are large intervals in 
which g is constant, intercalated by ranges of b where g varies wildly. If 6o li^s within one of 
these latter ranges, the final outcome of the light ray is highly uncertain. In Fig. 4b we show 
a magnification of one of these regions. Except for the scale, it is very similar to Fig. 4a. 
A further magnification is shown in Fig. 4c, again revealing structure in small scales. We 
have obtained even further magnifications, which are not shown here, and all show similar 
structures, down to the smallest scales allowed by the numerical limitations. This shows 
that g has a fractal dependency on b. Notice that there are large intervals of b where g 
is perfectly regular. These regular regions are mixed in all scales with the fractal regions, 
where the outcome of a light ray is highly uncertain. This sensitivity of the dynamics to the 
initial conditions is made precise with the definition of the box-counting dimension, which 
we now present briefly [Q. 

We define the basin boundary of the system to be the set of points (initial conditions) 

9 



such that all neighborhoods of these points have points belonging to at least two different 
basins, no matter how small that neighborhood is. The fractal nature of the basins shown 
in Fig. 4 results from a fractal basin boundary [Q. It is not difficult to see that a fractal 
basin boundary implies a fundamental uncertainty in the final outcome of an orbit. We 
now define the box-counting dimension of the basin boundary, which gives a measure of 
this uncertainty. Let bo be a randomly chosen impact parameter in the interval [be, besc]', we 
consider 0o = throughout for simplification. Let /(e) be the probability that there is a 
point of the basin boundary lying within a distance e from bo. In the limit e — > 0, / generally 
scales with e by a power law. We thus write: 

/(e) ex e'-'. (10) 

d is the box-counting dimension of the intersection of the basin boundary with the one- 
dimensional section of initial conditions given by 6 G [be, besc] and 0o = 0. Clearly, we must 
have < d < 1. If the basin boundary is regular, then d = 0; fractal boundaries have d > 0. 
f can be interpreted as a measure of the uncertainty as to which basin the point b belongs, 
for a given error e in the initial condition, which is always present in a real situation. For a 
regular basin boundary, / decreases linearly with e. If we have a fractal boundary, however, 
the power in ( [lO|) is less than 1, and / decreases much more slowly with e, which makes the 
uncertainty in the outcome much higher than in the case of a regular boundary. Thus, d is 
a good measure of the sensitivity to the initial conditions that results from a fractal basin 
boundary, and since it is a topological invariant 0, it is a meaningful characterization of 
chaos in general relativity. 

We calculate the box-counting dimension d numerically by using the method we now 
explain [Q. We pick a large number of initial conditions b randomly, and for each one of 
them we compute the map (§-|^), finding out its outcome and therefore to which basin it 
belongs. We then do the same thing to the two neighboring initial conditions b + e and 
b — e, for a given (small) e, for each b. If the three points do not belong to the same basin, 
b is labeled an "uncertain" initial condition. For a large number of initial conditions, we 
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expect that the fraction of uncertain points for a given e approximates /(e). Calculating in 
this way / for several values of e, we use Eq. (^0]) to obtain d from the inclination of the 
log-log plot of / versus e. Applied to the two-black-hole map with D = 15M, this method 
gives d = 0.17 ± 0.02. The error comes from the statistical uncertainty which results from 
the finite number of points used in the computation of /. In our calculation, the number of 
initial conditions was such that the number of "uncertain points" is always higher than 200. 

How does the fractal basin boundary arise from the dynamics of the map (H-^? In order 
to answer this question, we first observe that every point in the basin boundary gives rise to 
orbits that neither escape nor fall into one of the black holes (otherwise they would be part 
of one of the basins, which violates the definition of the basin boundary); that is, the basin 
boundary is made up of "eternal orbits" which move forever around the two black holes. We 
need thus to investigate these orbits to understand the formation of the basin boundary. 

Consider the one-dimensional set of initial conditions parameterized by the impact pa- 
rameter b with 00 = 0. We have seen that if |6| < be the orbit falls into the event horizon 
of a black hole, and if |6| > besc the orbit escapes. Thus, the points of the basin boundary 
belong to the interval 

\b\E[bc,besc], (11) 

which is actually two disjoint intervals, corresponding to positive and negative values of b. 
However, not all points in this interval are part of the basin boundary, of course; in order 
to survive the next iteration of the map (0-^) without escaping or falling, the corresponding 
orbits must be deflected in such a way that they reach the other black hole with an impact 
parameter within the interval (|ll]). From Fig. 5 we see that for this to happen the orbits 
must be deflected by an angle 6 in the neighborhood of (2n + l)7r + a, and either {2n + l)7r 
or (2n + 1)77 + 2a, depending on the previous deflection suffered by the orbit; the angle a 
depends on the distance separating the black holes, n is the number of turns the orbit makes 
around one of the black holes before moving on to the other one. For each n, there are two 
intervals of the deflection angle 6 for which the orbit survives the next iteration without 
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escaping or falling; these two intervals correspond to the positive and negative values of b 



satisfying Eq. ([Tl|). In the first iteration, the initial interval (|Tl]) is divided into infinitely 
many pairs of sub-intervals, each pair labeled by the number n of times the orbit circles the 
black hole. From Eq. (|^, intervals corresponding to large n's decrease exponentially with 
n. In the next iteration, each of these sub-intervals are themselves divided into an infinite 
number of intervals, and so on in the next iterations. In the limit of infinite iterations, the 
set of surviving orbits is a fractal set with zero measure. This set is the basin boundary, 
and its fractality is responsible for the complex dynamics shown in Fig. 4. The two fractal 
regions on the right of Fig. 4a consist of orbits whose first scattering has n = 0, that is, 
they are deflected by the black hole by vr and tt + a. The leftmost fractal region in Fig. 4a 
is actually an infinite number of very small regions, corresponding to orbits with n ^ 0; the 
scale of Fig. 4a is too large for them to be distinguished. This gives us an indication that 
the orbits with n > are a very small fraction of the basin boundary; we shall return to 
this later in this Section. 

Each orbit which is part of the basin boundary can be labeled by an infinite sequence of 
symbols 010203 ■ ■ ■ (including bi-infinite sequences • ■ ■ a_2a_iaoaia2 ' ' '), where each symbol 
Om = nm{km) givcs the neighborhood of the deflection angle {2nm + l)vr + fcm« after the 
m-th scattering. As an example, the periodic orbit shown in Fig. 3a is represented by 
the sequence 0(0)0(0)0(0) ■ ■ ■. In general, periodic orbits correspond to repeating sequences. 
However, not all sequences are allowed. It is clear from Fig. 5 that a symbol n(0) must be 
followed either by one of type m{0) or m{l), but it cannot be followed by a symbol like m{2), 
that is, of the form m{k) with k = 2. Analogously, a symbol with k = 1 oi k = 2 cannot be 
followed by one with k = 0. Even with these restrictions, however, there is an uncountable 
set of non-repeating sequences which label orbits that are part of the basin boundary. The 
uncountability of this set is a reflection of the fractal nature of the boundary. 

The basin boundary is the stable manifold of the invariant set, which is made by orbits 
labeled by bi-infinite symbols ■ ■ ■a_2a-iaoaia2 ■ ■ ■. These are orbits that do not escape for 
both forward and backward iterations of the map (H-^. 

12 



It is important to observe that the basin boundary is fractal because the scattering 
function A(6) of the isolated black hole (^ assumes values higher than tt, which makes 
it possible for orbits to be scattered to both sides of the black hole, giving rise to the 
fractal basin boundary. The scattering of particles by two fixed Newtonian mass points is 
immediately seen to be regular, because Rutherford's scattering function does not assume 
values higher than vr. This is of course in accordance with the fact that the fixed two- mass 
problem in Newtonian gravitation is integrable, since the Hamilton- Jacobi equation of this 
system is separated in elliptical coordinates |P|. 

We have seen that after being scattered by a black hole, the light ray must have an 
impact parameter lying on the interval (|TT]) to belong to the basin boundary. Since be = 
3\/3M ^ 5.19615M and besc ~ 5.35696M, the impact parameter must belong to one of two 
intervals of length Ah = b^sc — be ^ 0.16081M (the two intervals correspond to positive and 
negative impact parameters). In our approximation, we have Ab ^ D. Using this fact, we 
can approximate the value of b in Fig. 5 by besc, with an error of Ab at worst, which means 
a fractional error of about Ab/besc ~ 0.03. The distance L in Fig. 5 traveled by orbits which 
were deflected by (2n + l)7r + a or {2n + l)7r + 2a in the previous scattering is thus given 
by (L/2)2 + ble = {D/2)\ that is. 



D^ - 46L- (12) 

For the map (0-^) to be well-defined, we must have D > 2besc- Of course, this condition 
is satisfied in our approximation D ^ 2M. The angle a is calculated from Fig. 5 in this 
approximation, using elementary geometry: 



sma 



(13) 



D 

Consider a set of orbits with impact parameters filling the interval [6c, &esc] of length Ab, 
which may have already suffered several previous scatterings. The subsets of these orbits 
that survive the next scattering without escaping nor falling in one of the black holes are 
sub- intervals in the neighborhood of 6j^, where 6^ are the values of the impact parameter 
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such that the orbit is deflected by an angle of (2n + l)7r + ka {k = 0, 1 or 2). They are 
solutions of the algebraic equation 

A(6^) = (2n + l)7r + ka, (14) 

with A given by Eq. @. The fractal regions of Fig. 4a are located around values of b given 
by Eq. (0) with k = and k = 1. Depending on the previous deflections suffered by a 
given set of trajectories, the values k = 1 and k = 2 must be used instead in Eq. (^4|), 
according to the rules of the symbolic dynamics we exposed above. 

Because the distance D is much larger than Ab, the allowed range in the deflection angle 
69'^ around (2n + l)7r + ka of an orbit such that it arrives at the other black hole with b in 
the interval (|1T]) is approximately: 

K = ^-, M];- = ^, (15) 

where L is given by Eq. (0). The length Afe^ of the interval of surviving orbits around 6^ is 
in this approximation much smaller than A6: A6,^ ^ A6. We can therefore approximate A6^ 
by the first-order expression A6^ ^ 56^/\/S.'{b^)\, with A' = dA/db. We define A^ = A6,^/A6 
as the fraction of the interval [6c, bf.sc\ (or [— fogso ~^c\) occupied by the surviving orbits with 
b around 6^. We have: 

"^ I)|A'(60)|' ^" L|A'(6i'')r ^ ^ 

Given a value of D, it is easy to solve Eq. (|T^) numerically for b\, with L and a given 
by Eqs. (jg) and (0). For D = 15M, we find L ^ 10.5M and a ^ 0.7956. By a direct 
numerical calculation of A and its first derivatives, we find: 

6° = 6e.c; A'(60) = -5.863/M; A"(6|]) = SS.SO/M^; (17) 

bl = 5.26629M; A'(6^) = -13.35/Af; A"(6j) = 202.5/M2; 

bl = 5.22729M; A'(62) = -31.66/M; A"(62) = lOSO/M^; 
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6? = 5.19643M; A'(6?) ^ -3600/M; A"(6?) ^ 1.3 x W/M"^. 

Here A" = dP^/dh^. The second derivative of A will be used below. Notice that 6° is 
independent oiD. Using Eq. (p!6|), we obtain: 



\l = 0.01137; Aj = 0.007134; A^ = 0.003008; 

A° ^ 1.85 X 10"^ 

We see that A? is two or three orders of magnitude smaller that Ag, and the other A^'s 
with n > 1 are even smaller. From Eq. (^, they decrease exponentially with increasing 
n. These values show that for most purposes we can disregard the contributions to the 
dynamics from deflections with n > 0: the measure of the set of orbits that make multiple 
turns around a black hole is negligible. This approximation will be used extensively in the 
next Section. 

Equations (|I4|-^) are not exact, because the derivative A'(6) varies in the intervals A6^. 
To estimate the error, we use the fact that the A6^ are small. The error 6\^ in A^ is then 
to first order: 



Using A6^ = A6A,^, we get the fractional error 5A,^/A^: 

SX: A6 A"(60) _ 5X1^' Ab A"{h];^) 



1,2\ 



2- 



A'{h'n''f 



(18) 



A" D [A'(60)]^' A^'^ 
The terms A"(6^)/ ( A'(6^)) are of the order of 1, and the terms Ab/D and Ah/ L are much 



smaller than 1. Thus, the fractional errors in the values A,^ given by the approximate formula 
(p!6|) are very small, of about 0.01 for D = 15M. In the limit of large D, we have L ^ D, 
and (5A^/A^ ~ 1/D: the fractional error is inversely proportional to the separation D, and 



the approximation (|T^) gets better and better as D increases. 
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IV. THE LIMIT D ^oo 

We now take the limit D ^ 2besc- From Eqs. (0) and (|l2l), we have that in this hmit 
a^OandL^ D. Eqs. (|14D and (|16D then imply that 6^ ^ 6° = 6„ and A6^ -^ A6° = A6„: 
the two intervals of surviving orbits of a given n come closer and closer as D increases, and 
their lengths become the same in this limit. Because of the approximate equality in the 
lengths, the magnification of each interval by (A„)^^ gives approximately the same set of 
intervals. In other words, the fractal basin boundary is self-similar in this approximation. 

The box-counting dimension of this self-similar set is given by the solution of the tran- 
scendental equation g: 

oo 
n=0 

As we have seen in the previous Section, Ai is many orders of magnitude smaller than Aq. 
Therefore, it is a good approximation to neglect terms with n > in the above expression, 
and d can then be explicitly written as: 

In 2 

From Eq. (|I6D, we have (Aq)^^ = D\A'{besc)\, and we get: 

In 2 

d = ^-P^ t;, (20) 

lnL) + /3' ^ ^ 

with (3 = In |A'(&esc)|- For large D, the box-counting dimension decays logarithmically 

with the distance, which is a very slow decay law. Even for the large distances found 

in astronomy, d is still non-negligible. Since the box-counting dimension is linked to the 



scattering cross-section |jTO|, the slow decay of d with D shown in Eq. ( poD could have 



observational consequences. 

V. THE CASE OF DIFFERENT MASSES 

Now we consider the case of two black holes with different masses Ma and Mb. We 
restrict ourselves to the limit D ^ 26^^^, 26^^^, where 6"^^ and b^^^ are the impact parameters 
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corresponding to deflections of vr of tlie two black holes. Because the two masses are different, 
the ranges Aba and Abb of impact parameters for surviving orbits are different, and therefore 
the allowed range of scattering angles depends on which black hole the orbit is heading to. 
This means that the 'shrinking factors' Aq and Aq (given by Eq. ([T6|)) depend on the black 
hole. This spoils the property of self-similarity, which is a feature of the equal-masses case 
in the limit D -^ oo. However, since the orbits are scattered alternately by the two black 
holes, the square of the scattering map defined by the system is self-similar, in the limit 
D ^ oo. 

After two iterations, of each interval [be, besc] there remains four subsets of surviving 
orbits, all with size of approximately AqAq (we are not considering orbits with n > 0). 
Each of these subsets gives rise to four others after two further iterations, and so on. The 
box-counting dimension of the surviving set is given by [§ A{XqXqY = 1, that is: 

In 4 

a -- 



In (AgA^) ■ 

Substituting Eqs. (|16|) and (p!7D with M replaced by Ma and Mb, we obtain Aq and Aq, and 
we find: 

In 2 
lnD + /3 + ln^' ^^^^ 

where rj = Ma/ Mb, (3 = lnA'(6^^J. The difference between Eqs. (|2T]) and (gUp is the 
constant term In y^ in the denominator of d. If Ma = Mb, then rj = 1 and Eq. (^Tj) is equal 
to Eq. (^), as of course it should. In the limit D ^ oo, d also decays logarithmically with 
distance in this case, as in the case of equal masses. 

VI. CONCLUSIONS 

In this article we have studied the chaotic behavior of light rays orbiting a system of 
two non-rotating fixed black holes. We have assumed that the black holes are sufficiently 
far away from each other, so that we could consider the motion of the light rays to be 
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the result of the action of each black hole separately. Since the equations of motion of a 
light ray in the space-time of an isolated black hole can be solved analytically, using this 
approximation we reduce the motion of the massless test particle to a 2-dimensional map. 
Numerical integration of this map showed the existence of a fractal basin boundary, with 
an associated fractional box-counting dimension. In the limit of a large separation distance 
D between the two black holes, we have been able to obtain an analytical expression to the 
asymptotic value of the box-counting dimension d. We found that d ~ (InD)^^ for large D; 
this result also holds for different black hole masses. 
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FIGURE CAPTIONS 

Figure 1 The deflection function A (6) for an isolated Schwarzschild black hole. 

Figure 2 Construction of the scattering map. 

Figure 3 Two examples of periodic orbits of the scattering map. 

Figure 4 Portrait of the basins as a function of the impact parameter b, for 0o = 0. The 
'basin function' g{b) is defined to be 1 if the orbit with initial conditions (po = 0, bo = b 
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falls into one of the black holes, -1 if it falls into the other black hole, and if it 
escapes to the asymptotically plane region. Even though the function g{b) assumes 
only discrete values, the points are connected by straight lines for better visualisation. 
Two successive magnifications of a small region of the initial interval are pictured, 
showing the fractal dependency of g on b. 

Figure 5 (a) Two possible types of scatterings for an orbit which was previously scattered 
by (2n+ l)7r; (b) two scatterings for an orbit which was previously deflected by (2n + 
1)71 + a or (2n+ l)7r + 2a. 
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